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Abstract 

We present a new class of integrators for stiff PDEs. These integrators are 
generalizations of FLow Averaging integratORS (FLAVORS) for stiff ODEs and 
SDEs introduced in ^2] with the following properties: (i) Multiscale: they are based 
on flow averaging and have a computational cost determined by mesoscopic steps 
in space and time instead of microscopic steps in space and time; (ii) Versatile: 
the method is based on averaging the flows of the given PDEs (which may have 
hidden slow and fast processes). This bypasses the need for identifying explicitly 
(or numerically) the slow variables or reduced effective PDEs; (iii) Nonintrusive: A 
pre-existing numerical scheme resolving the microscopic time scale can be used as 
a black box and easily turned into one of the integrators in this paper by turning 
the large coefficients on over a microscopic timescale and off during a mesoscopic 
timescale; (iv) Convergent over two scales: strongly over slow processes and in the 
sense of measures over fast ones; (v) Structure-preserving: for stiff Hamiltonian 
PDEs (possibly on manifolds), they can be made to be multi-symplectic, symmetry- 
preserving (symmetries are group actions that leave the system invariant) in all 
variables and variational. 



1 Introduction 

Multi-scale PDEs can be divided into two (possibly over-lapping) categories: PDEs 
with highly oscillating or rough coefficients and PDEs with large (or stiff) coefficients. 
Classical numerical methods are usually: (i) stable but arbitrarily inaccurate for the 
former category (consider, for instance, a finite element method for the elliptic oper- 
ator — div(aV) with a rapidly changing coefficient a G L°°), or (ii) unstable for the 
latter category. Accurate numerical methods for the former category, called numerical 
homogenization methods, are, in absence of local ergodicity or scale separation, based 
on the compactness of the solution space (we refer, for instance, to [26^ [2l [27]). Numer- 
ical methods for the latter category are, in essence, based on the existence of slow and 
fast variables (or components) [14J. When fast variables converge toward Dirac (single 
point) distributions, asymptotic-preserving schemes jlS] allow for simulations with large 
time steps. We also refer to [THl [23] for multi-scale transport equations and hyperbolic 
systems of conservation laws with stiff diffusive relaxation. Well-identified slow variables 



can be simulated with large time-steps using the two-scale structure of the original stiff 
PDEs (we refer to pj and [jL2j for existing examples; slow variables satisfy a non-stiff 
PDE that can be identified in analogy to equations (A. 9) and (A. 13) of [32j; we also 
refer to [13] for a definition of slow variables). 

In this paper, we consider the second category of PDEs and propose a generaliza- 
tion of FLow Averaging integratORS (FLAVORS) (introduced in [32j for stiff ODEs 
and SDEs) to stiff PDEs. Multi-scale integrators for stiff PDEs are obtained without 
the identification of slow variables by turning on and off stiff coefficients in single-step 
(legacy) integrators (used as black boxes) and alternating microscopic and mesoscopic 



time steps (Subsection 2.2). We illustrate the generality of the proposed strategy by ap- 
plying it to finite difference methods in Section [2| multi-symplectic integrators in Section 
[sj and pseudospectral methods in Section [4] (although we have not done so in this paper, 
the proposed strategy can also be applied to finite element methods or finite volume 
methods). The convergence of the proposed strategy, after semi-discretization in space. 



is analyzed in Subsection 5.1, where a non-asymptotic error bound indicates the two- 
scale convergence ([32], i.e., strong with respect to hidden slow variables and weak with 
respect to hidden fast variables) of PDE-FLAVORS. As illustrated by numerical (Figure 
2]) and theoretical results (Section [5]), an explicit tuning ((/i/e)^ <C -ff <C h/e) between 
microscopic h and mesoscopic {H) time-steps and the stiff parameter 1/e is necessary 
and sufficient for convergence. We also show in Section [6] that applying the FLAVOR 
strategy to characteristics leads to accurate approximations of solutions of stiff PDEs. 

These results, along with those of [32j, diverge from the concept that, in situations 
where the slow variables are not linear functions of the original variables, multiscale 
algorithms "do not work" "if the slow variables are not explicitly identified and made 
use of" (page 2 of [T3]). 



2 Finite difference and space-time FLAVOR mesh 

2.1 Single-scale method and limitation 

Consider a multiscale PDE: 

F{1, X, t, u{x, t),Ux{x, t),ut{x, t),Uxx{x, t),Uxt{x, t),utt{x, t), . . .) = (1) 

where F is a given function (possibly nonlinear) , e is a small positive real parameter and 
X and t are spatial and temporal coordinates. 

To obtain a numerical solution of Q, the simplest single-scale finite difference ap- 
proach employs a uniform rectangular mesh with time step length h and space step length 
k, and approximates the solution u at its values at discrete grid points. Differential oper- 
ators will be approximated by finite differences; for instance, according to forward space 
forward time rules: Ux{ik,jh) ~ (tij+ij —Uij)/k and ut{ik,jh) ~ ("Uij+i —Uij)/h, where 
Uij is the numerical solution at discrete grid point with space index i and time index 
j. After this discretization, the original PDE is approximated by a finite dimensional 
algebraic system, which can be solved to yield the numerical solution. 
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Of course, a necessary condition for obtaining stability and accuracy in the numerical 
solution is that h and k have to be small enough. A quantitative statement on how 
small they need to be will depend on the specific PDE and discretization. For ID 
linear advection equations Ux — aut = and forward time forward space discretizations, 
the h < k/a CFL condition has to be met to ensure stability, which is also a 
neccessary condition for accuracy [TH]. Intuitively, the CFL condition guarantees that 
information does not propagate faster than what the numerical integrator can handle. 
The Von Neumann stability analysis [9j helps determine analogous CFL conditions for 
linear equations with arbitrary discretizations. The stability of numerical schemes for 
general nonlinear equations remains a topic of study. We refer to [31] for additional 
discussions on single-scale finite difference schemes. In general, the presence of a stiff 
coefficient e^^ in equation ([T]) requires h and k to scale with e in order to guarantee the 
stability of numerical integration schemes. This makes the numerical approximation of 
the solution of ([T]) computationally untractable when e is close to 0. 

2.2 Multiscale FLAVORization and general methodology 

FLAVORS are multiscale in the sense that they accelerate computation by adopting 
both larger time and space steps. A finite difference scheme can be FLAVORized by 
employing two rules: 
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Fi gur6 1* Mesh used by FLAVORS. A uniform mesoscopic space step is used and two alternating small 
and mesoscopic time steps are used. Stiffness is turned on in red regions and turned off otherwise. 
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First, instead of a uniform mesh, use a mesh as depicted in Figure [T| in which a 
uniform spatial grid corresponds to a mesoscopic space step K that does not scale with 
e, and an alternating temporal grid corresponds to two time steps, microscopic h (scaling 
with e) and mesoscopic H — h (H independent from e). It is worth mentioning that when 
using this non- uniform mesh, grid sizes have to be taken into consideration when deriva- 
tives are approximated by finite differences. Ist-order derivatives are straightforward to 
obtain, and we refer to Section [3] for approximations of higher order derivatives. 

Second, the stiff parameter should be temporarily set to be (i.e., turned off) 
when the current time step is the mesoscopic H — h; if the small time step h is used 
instead, the large value of needs to be restored, or in other words, stiffness should 
be turned on again. 

The rule of thumb is that k and h should be chosen such that the integration of ([T]) 
with these step sizes and stiffness turned on is stable and accurate. On the other hand, 
there is another pair of step size values such that the same integration with stiffness 
turned off is stable and accurate, and K and H should be chosen to be an order of 
magnitude smaller than these values. FLAVORS does not require a microscopic k, but 
only a mesoscopic space-step K, a microscopic time-step h, and a mesoscopic time-step 
H. 

The intuition is as follows: adopt the point of view of semi-discrete approach for 
PDE integration, in which space is discretized first and the PDF is approximated by a 
system of ODFs. The integration (in the time) of the resulting finite dimensional ODE 
system can be accelerated by applying the FLAVOR strategy to any legacy scheme 
(used as a black box). Turning on and off stiff coefficients in the legacy scheme and 
alternating microscopic time steps (stiffness on) with mesoscopic time steps (stiffness 
on) preserves the symmetries of that scheme and at the same time induces an averaging 
of the dynamic of (possibly hidden) slow variables with respect to the fast ones. With 
this strategy, the FLAVORized scheme advances in mesoscopic time steps without losing 
stability. The (possibly hidden) slow dynamic is captured in a strong sense, while the fast 
one is captured only in the (weak) sense of measures. A rigorous proof of convergence 
of the proposed method relies on the assumption of existence of (possibly hidden) slow 
variables and of local ergodicity of (possibly hidden) fast variables (we refer to Section [s]). 
It is important to observe that the proposed method does not require the identification 
of slow variables. 

2.3 Example: conservation law with Ginzburg-Landau source 

Consider a specific stiff PDE: 

ut + f{u), = e-\il-u^) (2) 

in which f(u) = s'mu and < e ^ 1. Use the boundary condition of u{x = 0,t) = 
u{x = L,t) and the initial condition of u{x,t = 0) = sin(7rj;). This system contains two 
scales: the fast process corresponds to u quickly converging towards 1 or —1, and the 
slow process corresponds to the front (with steep gradients) separating u > from u < 
propagating at an 0(1) velocity. 
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We will FLAVORize the following Lax-Friedrichs finite difference scheme: 



, "2+1 J — 2 



(3) 



where ttjj = u^^i/f^j and tij^i = sin (7r(i — 1)A;). If the domain of integration is restricted 
to [0, L] X [0, T], then i = 1, 2, . . . , [L/A;J +1, and j = 1, 2, . . . , [T/h\ +1. We use h = O.le 
and = 0.2e for our purposes, both of which we found numerically at the order of the 
stability limit. In our experiment, we chose e = 2 • 10~^, and therefore h = 0.0002 and 
k = 0.0004. 

The FLAVORized version of this scheme is: 



Ui+ij = {ui+2,j + Uij)/2 (4) 



2 / 2K 



where Uij = "Ui+L/i^-j and Uj^i = sin (7r(i — 'i)K). If the domain of integration is re- 
stricted to [0, L] X [0,'t], then i = 1, 2, . . . , [^/i^J + 1, and j = 1, 2, . . . , [T/H\ + 1. We 
use the same h as before, and choose H = 0.005 and K = 0.01, which ensures that the 
stability of the integration remains independent of e. 

Errors of FLAVOR based on Lax-Friedrichs with different H and h values are com- 
puted by comparing the results to a benchmark Lax-Friedrichs integration with fine steps 
h = O.le and k = 0.2e. More precisely, we calculated the distance between two vectors 
respectively corresponding to FLAVOR and Lax-Friedrichs integrations, which contain 
ordered m(x, t) values on the intersection of FLAVOR and Lax-Friedrichs meshes (which 
is in fact the FLAVOR mesh as long as is a multiple of O.le). 1-norm is used and 
normalized by the number of discrete points to mimic the norm for the continuous 
solution. Experimental settings are e = 2 • 10~^ L = 2 and T = 2. As we can see 
in Figure [2j FLAVOR is indeed uniformly convergent in the sense that the error scales 
with as long as h takes an appropriate value. This is not surprising, because we have 
already proven in the ODE case that the error is bounded by a function of H (uniformly 
in e) as long as (^) <^ H <^ h/ e, and this error can be made arbitrarily small as H 
(notice H can still be much larger than e as e J, 0). 

Also, a typical run of FLAVOR {H = 0.005 and K = 0.01) in comparison to the 
benchmark {h = 0.0002 and k = 0.0004) is shown in Figure [sj FLAVOR captured 
the slow process strongly in the sense that it obtained the correct speeds of both steep 
gradients' propagations (up to arithmetic error and fringing). In this setting, FLAVOR 
achieves a = 312.5 fold acceleration. It is worth restating that both spatial and 
temporal step lengths of FLAVOR are mesocopic, whereas the counterparts in a single 
scale finite difference method have to be both microscopic for stability. The computa- 
tional gain by FLAVOR will go to infinity as e — )• 0, and this statement will be true for 
all FLAVOR examples shown in this paper. 
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L, distance between u . by FLAVOR and benchmark 
1 x,t ' 





Figure 2: Errors of FLAVOR based on Lax-Friedrichs as a function of H and h. H samples multiples 
of O.le, starting from 2x to 50x with Ix increment, and h ranges from O.Ole to 3e with O.Ole increment. 
Errors with magnitude bigger than 1 are not plotted, for they indicate unstable integrations. 



3 Multisymplectic integrator for Hamiltonian PDEs 

3.1 Single-scale method 

We refer to 13 [2T| [22] for a discussion on the geometry of Hamiltonian PDEs (e.g., 
multi-symplectic structure). We will now recall the Euclidean coordinate form of a 
Hamiltonian PDE: 

Mzt + Kz^ = V,H{z) (5) 

where z(x, t) is a n-dimensional vector, Ai and fC are arbitrary skew-symmetric matrices 
on M", and H : ^ M is an arbitrary smooth function. The solution preserves the 
multi-symplectic structure in the following sense: 

dti{u, V) + d^u, v) = o (6) 

where l and k are differential 2-forms defined by 

i{x,y) = {Mx,y) and K{x,y) = {ICx,y) (7) 
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Single scale finite difference FLAVOR 




and U and V are two arbitrary solutions to the variational equation (the solution is 
identified with dz : R"^ ^ W): 

Mdzt + fCdz^ = D^^H{z)dz, dz{x, t) e (8) 

Preservation of multi-symplecticity can be partially and intuitively interpreted as a con- 
servation of infinitesimal volume in the jet bundle, which generalizes the conservation of 
phase space volume in Hamiltonian ODE settings to field theories. 

A broad spectrum of PDEs fall in the class of Hamiltonian PDEs, including general- 
ized KdV, nonlinear Schrodinger models, nonlinear wave equations, atmospheric flows, 
fluid-structure interactions, etc. jUISl El [7| . We also refer to [8j and references therein 
for surveys on numerical recipes, and to [20] for an application to numerical nonlinear 
elastodynamics . 

Hamiltonian PDEs ^ can be viewed as Euler-Lagrange equations for field theories, 
which are obtained by applying Hamilton's principle (i.e., a variational principle of 
6S/6z = 0) to the following action: 

S{z{-,-)) = jj C{z,zt,z^)dtdx (9) 

where the Lagrangian density is given by 

£{z, zt, z^) = ^{Mzt, z) + ^(/Cz^, z) - H{z) (10) 

This variational view of Hamiltonian PDEs will intrinsically guarantee the preserva- 
tion of multi-symplecticity, and there will be a field generalization of Noether's theorem, 
which ensures conservation of momentum maps corresponding to symmetries. 

Numerically, instead of discretizing the equations ([s]) , we prefer the approach of varia- 
tional integrators because they are intrinsically multi-symplectic and therefore structure- 
preserving |2 H 1221 [23l I20j. These integrators are obtained as follows: first discretize the 
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action ^ using quadratures, then apply variational principle to the discrete action 
(which depends on finitely many arguments), and finally, solve the algebraic system 
obtained from the variational principle, i.e., the discrete Euler-Lagrange equations. 
For an illustration, consider a nonlinear wave equation: 



Utt - Uxx 



V'(u) 



ill) 



with periodic boundary condition u{x + L,t) = u{x, t) and compatible initial conditions 
u{x, t = 0) = f{x) and ut{x, t = 0) = g{x). Suppose we are interested in the solution in 
a domain [0,L] x [0,r]. 

Rewrite the high order PDE as a system of first order PDEs (notice these covariant 
equations can be obtained through an intrinsic procedure, which works on manifolds as 
well [5]): 



Vt - Wx 
Ut 
Ux 



V 

w 



The corresponding Lagrangian density is: 

C 



1 2 1 2 

2 i 2 ^ 



V{u) 



(12) 
(13) 
(14) 



(15) 



Using a forward time forward space approximation, we obtain the following discrete 
Lagrangian: 



L 



h'lj h'lj 



1 f UiJ+l - Ui^j 

2 



hij 



1 f Ui+ij - u, 



tj+i-tj+hij 



dt 



dx 



1 



-u. 
2 * 



ul + V[u) 



(16) 
(17) 



where space step kij and time step hij define a rectangular grid of size kij x hij. The 
simplest single-scale choice would be kij = k and hij = h for some k and h. 

As a consequence, the continuous action S is approximated by a discrete action: 



N M 
a=l 13=1 



and Hamilton's principle of least action dSd = gives 

N M 



C dt dx 



(18) 



du 







o=l 13=1 



for 1 < i < and 1 < j < M, where N and M are such that XloLi 
and YlsLi ha (3 = T for any a. 



(19) 

L for any /? 
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Taking derivative with respect to Uij, we obtain the following discrete Euler-Lagrange 
equations: 



k 



-+hijkijV'{uij)+ki 



n-i-ij — 



(20) 

The system of above equations is explicitly solvable when equipped with boundary 
conditions and initial conditions; for instance, below is a consistent discretization of the 
continuous version: 



Ui,2 = + hug {J2a=l 



(21) 



This numerical receipt is convergent. In fact, multi-symplectic integrators obtained 
from variational principles can be viewed as special members of finite difference methods, 
whose error analysis is classical. 

We wish to point out that the above procedure works for any Hamiltonian PDEs 
of form ([5]). Also, notice that high-order derivatives are dealt with in an intrinsic way 
regardless of whether mesh is uniform. 



3.2 FLAVORization of multi-symplectic integrators 

Now consider a multiscale Hamiltonian PDE 

M{1, e~^)zt + /C(l, e-^)z^ = V,H{1, e-'\z) 



(22) 



Any single-scale multi-symplectic integrator can be FLAVORized (to achieve com- 
putational acceleration) by using the following strategy: (i) Use the two-scale mesh 
illustrated in Figure [l| and (ii) Turn off large coefficients when taking mesoscopic time- 
steps. Unlike FLAVORizing a general finite diff'erence scheme, we FLAVORize the action 
Sd instead of the PDE. Specifically, choose 



ki 



hi 



K, 



Vi, j 

Vi and odd j 



(23) 



H — h, Vi and even j 



in Lfj for even j's and all i's, while the large value of e ^ is kept in Lfj 



and let 

for odd j's and all z's. h and H correspond to a small and a mesoscopic time-step, and 
K corresponds to a mesoscopic space-step; the same rule of thumb for choosing them in 
Section [2] applies. 

After applying the discrete Hamilton's principle, the resulting discrete Euler Lagrange- 



equations corresponding to a multi-symplectic integrator will still be (20), except that 



stiffness is turned off in half of the grids. Multisymplecticity is automatically gained, 
because the updating equations originate from a discrete variational principle [21j. 
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3.3 Example: multiscale Sine-Gordon wave equation 



Consider a specific nonlinear wave equation (11) in which V{u) = — cos(u;u) — cos(ti). 
If w = 0, this corresponds to Sine-Gordon equation, which has been studied extensively 
due to its soliton solutions and its relationships with quantum physics (for instance, as 
a nonlinear version of Klein-Gordon equation). We are interested in the case in which 
oj (identified with e~^) is big, so that a separation of timescale exhibits. 

Arbitrarily choose L = 2 and use periodic boundary condition u{x + L,t) = u{x,t), 
and let initial condition be u{x,0) = sin(27rx/L) and ut{x,0) = 0. Denote total simula- 



tion time by T. Use the FLAVOR mesh (23). In order to obtain a stable and accurate 



numerical solution, k and h have to be o{l/uj), and K and H need to be o(l). 



Single scale 1st order multlsymplectic: 01=20 



FLAVOR:o)=20 













Figure 4: Numerical solutions to multiscale Sine-Gordon equation by single-scale Ist-order multi- 
symplectic integrator (left) and its FLAVORization (right). For clarity, the surface plots (but not 
simulations) use the same mesh size. 

A comparison between the benchmark of the single-scale forward time forward space 
multi-symplectic integrator (Eq. l20]with hij = h and kij = k) and its FLAVORization 



(Eq. [20] with mesh (23) and V'{u) = ujs'm{uju) -|- sin(n) for odd j and V'{u) = sin(n) 
for even j) is presented in Figure [4j uj = 20, k = L/20/lo and h = k/2, and K = L/40 
and H = K/2. It is intuitive to say that the slow process of wave propagation is well- 
approximated by FLAVOR, although the fast process of local fluctuation is not captured 
in the strong sense. Error quantification is not done, because what the slow and fast 
processes are is not rigorously known here. HK/2hk = 50-fold acceleration is obtained 
by FLAVOR. 

Readers familiar with the splitting theory of ODEs [24J might question whether 
FLAVORS are equivalent to an averaged stiffness of u) = uij^ (which corresponds oj = 2 
in the numerical experiment described above). The answer is no, because the equivalency 
given by the splitting theory is only local. In fact, the same single-scale forward time 
forward space multi-symplectic integration of the case a; = 2 is shown in Figure [5| which 
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Single scale 1st order multisymplectic: (o-2 




Figure 5: Numerical solutions to multiscale Sine-Gordon equation with the 'equivalent' stifTness by 
single-scale Ist-order multi-symplectic integrator. For clarity, the surface plot (but not the simulation) 
uses the same mesh size (as in Figure |4|. 

is clearly distinct from the FLAVOR result in Figure |4| Moreover, because of the e'^^'^ 
error term, changing stiffness alone will not result in a converging method (and result 
in a 0(1) error on slow variables). 



4 Pseudospectral methods 

4.1 Single-scale method 

Consider a PDE 

ut{x,t) = Cu{x,t) (24) 

with periodic boundary condition u{x, t) = u{x + L, t) and initial condition u{x, 0) = 
/(x), where i2 is a differential operator involving only spatial derivatives. 

The Fourier collocation method approximates the solutions by the truncated Fourier 
series: 

UM{x,t)= a„(t)e^"2Wi (25) 

|n|<Ar/2 

and solves for a„(t)'s by requiring the PDE to hold at collocation points yj: 

dtUN {Vj ,t) - Cum {yj,t) = (26) 

This yields N ODEs, which can be integrated by any favorite ODE solver. Of course, 
specific choices of collocations points will affect the numerical approximation. Often- 
times, the simplest choice of yj = Lj/N,j = 0, . . . , A'' — 1 is used, and in this case, the 
method is also called a pseudospectral method. We refer to [T7| for additional details on 
Fourier collocation methods. It is worth mentioning that pseudospectral methods can 
also be multi-symplectic when applied to Hamiltonian PDEs |10j . 
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4.2 FLAVORization of pseudospectral methods 

When the PDE is stiff (for instance, when L contains a large parameter e~^), FLAVORS 
can be employed to integrate the stiff ODEs (which will still contain e~^) resulting from 
a pseudospectral discretization. 

Similarly, for the FLAVORization of a pseudospectral method, it is sufficient to 
choose N ^ L instead of ^> e~^L, i.e., the space-step can be coarse {K = o(l)). 
For time stepping, alternatively switching between h = o(e) and H — h for a mesoscopic 
H = 0(1) is again needed, and stiffness has to be turned off over the mesoscopic step of 
H — h. In a sense, we are still using the same FLAVOR 'mesh' (Figure [T]), except that 
here we do not discretize space, but instead truncate Fourier series to resolve the same 
spatial grid size. 

4.3 Example: a slow process driven by a non-Dirac fast process 

Consider the following system of PDEs 

Ut + Ux- =Q 

qt + qx-p = (27) 
,Pt +Px + i^'^q = 

with periodic boundary conditions u{x,t) = u{x+L,t), q{x,t) = q{x+L,t), and p{x,t) = 
p{x + and initial conditions u{x,Q) = /"(x), g(x,0) = fix), and p{x,0) = P{x). 
The integration domain is restricted to [0, T] x [0, L]. The stiffness is identified with 
w^. We choose the initial condition of /"(x) = f'^{x) = cos{2itx / L) and P\x) = 0. 

In this system, q and p correspond to a fast process, which is a field theory version 
of a harmonic oscillator with high frequency oj. u is a slow process, into which energy is 
pumped by the fast process in a non-trivial way. 



We have chosen to FLAVORize (27) because it does not fall into the (simpler) cat- 
egory of systems with fast processes converging towards Dirac (single point support) 
invariant distributions |15j . 

We use the classical 4th order Runga-Kutta scheme (see, for instance, |16J) for the 



(single-step) time integration of the pseudospectrally discretized system of ODEs (26) 



Write (j)^ : an'^'^{t) 1— )• dn^''^{t + h) its numerical flow over a microscopic time step 



h (consisting of four sub-steps), where a^''''^(t) are numerical approximations to the 
Fourier coefficients in ( |25| ), for the unknowns u, q and p at an arbitrary time t. Then, 
the corresponding FLAVOR update over a mesoscopic time step H will be 4'%_j^ o cj)^ , 
which consists of eight sub-steps. 

We present in Figure [6] and Figure [7] a comparison between the benchmark of single- 
scale pseudospectral simulation and its FLAVORization. It can be seen that the slow 
process of u is captured in strong (point-wise) sense, whereas the fast process of q is only 
approximated in a weak sense (i.e. as a measure, in the case, wave shape and amplitude 
are correct, but not the period). We choose L = 2, T = 10 and to = 1000. The single- 
step integration uses = 20 and h = 0.1/tJ (notice that this is already beyond the 
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Pseudospectral: slow process 



FLAVOR: slow process 





Fi gurc 61 Single-scale (left) and multiscale pseudospectral (right) integrations of slow u in system ( |27| 
Plotting mesh for the single-scale simulation is coarser than its computation mesh. 



stability/accuracy region of a single-scale finite difference, since the space step does not 
depend on l/cj; the spectral method is more stable/accurate for a large space-step), and 
FLAVOR uses iV = 20, /i = l/w^ and H = 0.01. H/2h = 50-fold acceleration is achieved 
by FLAVOR. 



5 Convergence analysis 
5.1 Semi-discrete system 

All FLow Averaging integratORS described in previous sections are illustrations of the 
following (semi-discrete) strategy: first, space is discretized or interpolated; next, spatial 
differential operators are approximated by algebraic functions of finitely many spatial 
variables; finally, the resulting system of ODEs is numerically integrated by a corre- 
sponding ODE-FLAVOR [32j . In this section, we will use the semi-discrete ODE system 
as an intermediate link to demonstrate that these PDE-FLAVORS are convergent to 
the exact PDF solution under reasonable assumptions (in a strong sense with respect 
to (possibly hidden) slow variables and in the sense of measures with respect to fast 
variables) . 

More precisely, consider a spatial mesh (vector) Jv[^ = [xi,X2, • • •], a temporal mesh 
(vector) M.'^ = [ti,t2) • • •]) and a domain mesh (matrix) M = M.^ x M.'^ . Examples of 
these meshes include the FLAVOR mesh = [K, 2K, NK] and M'^ = [h, H,H + 
h, 2H, . . . , (M — l)H, (M — l)H + h, MH], and a usual single-scale (step) integration 
mesh = [k, 2k,.. . ,L] and M"^ = [h, 2h,. . .,T] (recah the domain size is L = NK 
by T = MH). We will use the FLAVOR mesh throughout this section. We will compare 



the solution of the PDE (28) with the solution obtained with the FLAVOR strategy at 
these discrete points. 
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Pseudospectral: fast process 



FLAVOR: fast process 




Figure 7: Single-scale (left) and multiscale pseudospectral (right) integrations of fast q in system \27\ . 
Plotting mesh for the single-scale simulation is coarser than its computation mesh. The same color does 
not indicate the same value in these two plots. 



For simplicity, assume the PDE of interest is Ist-order in time derivative: 
ut{x,t) = F{l,€''^,x,t,u{x,t),Uxix,t), . . .) 



(28) 



Observe that a PDE ([T]) with higher-order time derivatives can be written as a system 
of Ist-order (in time derivatives) PDEs. 

Now consider a consistent discretization of PDE (28) with space step K and time 
step h (we refer to Page 20 of |31j for a definition of the notion of consistency, which 
intuitively means vanishing local truncation error). Letting /i J, in this discretization, 
we obtain a semi-discrete system (continuous in time and discrete in space). This semi- 
discrete system is denoted by the following system of ODEs, with approximated spatial 
derivatives: 

uiit) = /i(ni,n2, . . .,UN,e~^,t) 

U2{t) = /2(W1,U2, . . .,UN,e'^,t) ^^g^ 
UN{t) = fN{ui,U2, . ■ .,UN,e~^,t) 

Assuming existence and uniqueness of an exact strong solution u to the PDE (28), 
and writing u{A4f,t) its values at the spatial discretization points, we define for each i 
the following remainder: 

M^~\t) = |^(-^f ' i) - fM^f, t),u{Ml t),..-, u{M%,t),e-\t) (30) 

which is a real function of t indexed by e^^. 

Then, Ui{t) approximates the exact solution u{A4^ ,t) evaluated at grid points in the 
sense that these remainders vanish as e~^K \. (where K := A4f — A4f_i): 
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Lemma 5.1. Assume that F in (28) satisfies 

\F{l,e-\x,t,u{x,t),u-^{x,t),...)\ < {I + e~^)\F{l,l,x,t,u{x,t),u-^{x,t), . . .)\ (31) 



Assume that the fi in ([29j) satisfy similar inequalities. Then, there exists a constant Ci 
independent from e, h, H or K, such that for hounded t and u 



\ni{e-\t)\ < {l + e-^)CiK 



(32) 



Remark 5.1. (31) is true, for instance, in cases where 



F{l,e ^ ,x,t,u{x,t), . . .) = FQ{x,t,u{x,t), . . .) + e ^Fi{x,t,u{x,t), . . .). (33) 



Proof. The Hnear scaUng with K in (32) immediately follows from the definition of 
consistency, and the parameter 1 + in (32) has its origin (31). □ 

Remark 5.2. The consistency of finite difference methods can be easily shown using 
Taylor expansions. For instance, applying a Taylor expansion to the solution of ut — 
e'^Ux = a{u) leads to 

/-T^ .,x , / ^} /uUi + l)K,jh) — u(iK,jh) 

^ K (34) 



which implies 
d 



u{iK, t) 



+ 0{K)) + a{u{iK,jh))) + 0{h^ 
_^u{{i + l)K,t) - u{iK,t) 



+ a{u{iK, t)) + e-^0{K) 



dt^'^"'^' ^ K 
and naturally establishes the correspondence of fi{ui, . . . ,UN,e~^,t) = g-i "'+iW~'"'(*) 



(35) 



K 



+ 



a{ui{t)) and TZi = €~^0{K) for a Ist-order finite difference scheme. Notice that the 
remainders are still stiff, but we will see later that this is not a problem, since they 
can be handled by ODE-FLAVORs. The consistency of pseudospectral method can be 
shown similarly using Fourier analysis. 



With TZi defined in (30), consider the following system of ODEs: 



ui{t) = fi{ui,U2,...,UN,e ^^) + 7^l(e ^^) 



(36) 



^Ujy{t) = fN{ui,U2, ■ ■ . ,UN,e '^,t)+nN{e \ *) 



with initial condition nj(0) = u{M.^, 0). Obviously, its solution {ui{t))i<i<N is the exact 
PDF solution sampled at spatial grid points, i.e., Ui{t) = u{M.^ ,t). 

We will now establish the accuracy of PDE-FLAVOR by showing that an ODE- 
FLAVOR integration of (36) leads to an accurate approximation of ('Ui(t))i<j<Ar. Since 
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space (with fixed width L) is discretized by N grid points, we use the foUowing (nor- 
mahzed by A^) norm in our following discussion (suppose Vi{t) = v{A4f , t) for a function 
v): 

\\[vi{t),V2{t),...,VNm\ = ^\\[vi{t),V2{t),...,VN{m\, (37) 

Observe that if v{-,t) is Riemann integrable, then 

1 



hm II [v{Mf,t),v{M^,t), v{M%,t) 



L 



v{-,t)\\^i (recall L = iVET is fixed), 

(38) 



and hence the norm (37) does not blow up or vanish as N ^ oo. 



5.2 Sufficient conditions for convergence, ODE-FLAVORS, and two- 
scale convergence of PDE-FLAVORS 

We will now prove the accuracy of PDE-FLAVORs under the assumption of existence 
of (possibly hidden) slow and locally ergodic fast variables. The convergence of PDE- 
FLAVORs will be expressed using the notion of two-scale flow convergence introduced 
in [32] (corresponding to a strong convergence with respect to slow variables and weak 
one with respect to fast ones). 



Condition 5.1. Assume that the ODE system (36) satisfies the foUowing conditions: 

1. (Existence of hidden slow and fast variables): There exists a (possibly time- dependent) 
diffeomorphism rf : [ni(t), . . . , UN{t)\ i-)- [x{t),y{t)\ from onto M^^*' x W with 
uniformly bounded C"^ , derivatives with respect to Ui 's and t, and such that for 
all e > 0, {x{t),y{t)) satisfies 




f{x{t),y{t),t) 
e-^g[x{t),y{t),t) 



where f and g have bounded derivatives with respect to x, y and t. 

2. (Local ergodicity of vast variables): There exists a family of probability measures 
^*(x, dy) on MP indexed by x ^ and t and a family of positive functions 

T I— >• E^(T) satisfying liniT-^oo (T) = for all bounded t, such that for all 
xo,yQ,tQ,T bounded and (p uniformly bounded and Lipschitz, the solution to 

Yt = 9{xo,Yt,to) Yo = yo (40) 

satisfies 



^cl>iYs)ds- I 0(y)M*«(xo,dy) < x*°(||(^o,yo)||)^*°(T)(||0||L^ + ||V</.|h 



1 

f 

(41) 

where r i— )■ X^°i''") bounded on compact sets, and ^* has bounded derivative with 
respect to t in total variation norm. 
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Under Conditions 5.1 the computation of the solution of PDE (28) can be acceler- 
ated by applying the FLAVOR strategy to a single-scale time integration of its semi- 



discretization (29). 



Write the numerical flow of a given (legacy) ODE integrator for (29): 



t,t+T 



: [ui{t), UNit)] ^ [ui{t + t), . . . , UNit + t)] , 



(42) 



where Ui{s) approximates Ui{s) for all s, r is the integration time step, and a is a 
controllable parameter that replaces the stiff parameter in (29) and takes values of 
(stiffness 'on') or (stiffness 'off'). 

Definition 5.1 (ODE-FLAVORS). The FLow Averaging integratOR associated with 
$ is defined as the algorithm simulating the process: 

[ui{t),...,UN{t)] 



{k-l)H+h,kH ° ^{k-l)H,{k-l)H+h} 



(43) 



where (the number of steps) /c is a piece-wise constant function of t satisfying kH < 
t < (k + l)H, /i is a microscopic time step resolving the fast timescale {h <^ e), H is a 
mesoscopic time step independent of the fast timescale satisfying h e <^ H <^ 1 and 



/^n2 tt h 

{-y <^H 

e e 



(44) 



Condition 5.2. Consider the legacy ODE integrator with one-step update map ^j^^+t- 
introduced in (42). Suppose there exists constants C > and Hq > independent of N 
and a, such that for any r < Hq min(l/a, 1) and bounded vector [ui, . . . , un], 



\<^f^t+^iui,...,UN) - [ui,...,un] -T[fi{ui,...,UN,a,t),... 

...,fN{ui,...,UN,a,t)]\\ < Cr'^il + af , 



Condition 5.2 corresponds to the assumption that the integrator 



(45) 

, t+T is consistent 

for ([29]). ' ' 

Observe that we are integrating (29) but not (36), since the remainders T^i's are a- 
priori unknown unless the exact PDE solution is known. However, the following lemma 
implies that the FLAVORization of this integration is in fact convergent to the solution 



of (36), even though T^j's are possibly stiff. 



Lemma 5.2. Assume that ^f^_^_^, introduced in (42), satisfies Condition 



5.2 



Let h 



and H be the time steps used in the FLAVORization 5.1 If h <^ e, H <C h/e, and 
K = 0{H), then 



■ ,un) - [ui, . . .,un] - r[fi{ui, . . .,UN,a,t) +TZi{a,t), 
fN{ui, . . . ,UN,a,t) +TZN{a,t)]\\ < CT^{l + af 



(46) 



where t = h when a = e ^ and t = H — h when a = 0. 
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Proof. By Condition |5.2[ we have 



\^t,t+riui, ■ ■ ■,un) - [ui, ...,un]- r[/i(ni, . . . ,UN,a,t), . . . 

...,fNiui,...,UN,a,t)]\\ < Cr^l + af (47) 



for any r < min(l/a, 1)Hq. In addition, Lemma 5.1 gives a bound on the remainders: 
when a = e~^, there exists a constant C > independent of N and e^^, such that for 
all i, 

\Tni{e-'^,t)\ <TCKe-^ (48) 

Because we use t = hin this case and K <^ e^^r, the above is bounded by TC{Ce~^T)e~^ < 
Cr^(l + a)^ for some constants (7 <C 1 and C = CC. When a = on the other hand, 
there exists a constant C > such that for all i 

|r7^i(e~^^)| < rCK (49) 

Because K = 0{H) and we use t = H — h = 0{H) in this case, the above is bounded by 
tCCt < Ct'^{1 + a)^ for some constants C and we let C = CC. Notice that the value 
of K is fixed in both cases but r has different values: the flow map used in FLAVOR 
associated with a = is the one with mesoscopic step i.e., t = H — h; when 

a = on the other hand, the flow map is ^If^f^ and t = h. Finally, the triangle 
inequality gives 

||$°t_^^(Mi, . . . ,nAr) - [ui,...,un] - T[fi{ui, . . . ,UN,a,t) +ni{a,t), . . . 
. . .Jn{ui, . . .,UN,a,t) +TZN{a,t)]\\ < ||$°(_|.^(ui, . . .,un) - [ui, . . .,un]- 

1 ^ 

r[/i(ni, . . .,UN,a,t), . . . ,/Ar(ni, . . .,UN,a,t)]\\ + ^ X] k'^i(">*)l < 2Cr^(l + af , 

i=l 

(50) 

which finished the proof after absorbing the coefficient 2 into C. □ 

We also need the usual regularity and stability assumptions to prove the accuracy of 
FLAVORS for ([36]). 

Condition 5.3. Assume that 

1- fi, f2, ■ ■ ■ , fN are Lipschitz continuous. 

2. For all hounded initial condition [ni(0), . . . , un{{))] 's, the exact trajectories 
([wi(t), . . . , nAr(t)])o<t<T (i-e., solution to (p6|) ) are uniformly hounded in e. 



3. For all bounded initial condition [ui{0), . . . ,un{0)]'s, the numerical trajectories 
{[ui{t), . . . ,UN{t)])o<t<T (defined hy (43)J are uniformly bounded in e, < H < 
Ho, h < mm{Hoe,H). 
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The following theorem shows the two-scale flow convergence (strong on slow variables 
X and in the sense of measures on fast ones y, see [32]) of FLAVORs under the above 
conditions. 



5.2 



and 



Th eore m 5 .1. Consider FLAVOR trajectories in Definition \5.1\ Under Conditions 5.1 
5.3. there exist C > 0, C > and Hq > independent from e^^ and N, such 



that for K/C < H < Hq, h < Hqc and t > 0, 

\\x{t) - [v'TiMt], UNimi < Ce^\i{ui{0), un{0), e, H, h) 
and for all hounded and uniformly Lipschitz continuous test functions ip : 



(51) 



1 

At 



t+At 



ip{[ui{s),...,UNis)])ds- / ip{[r]^] \x{t),y))fi\x{t),dy) 

JRP 

< X2{uiiO),...,UNiO),e,H,h,At,t)i\\ip\\L^ + ||Vv9||l^; 



(52) 



where xi o.i^d X2 cl^^ bounded functions converging towards zero as e < H/{Cln^), 
f i 0, j^H iO and (f )2^ i (and At i for X2)- 

Recall notations: NK = L is the fixed spatial width, [r/*]^ and respectively 
denote the x (slow) component and the inverse of the diffeomorphism (defined in 



Condition 5.1), x{t) = [r] ]^{ui{t), . . . ,MAr(t)) corresponds to the slow component of the 



exact PDE solution sampled at grid points. Ui{t) and Ui{t) represent the exact and the 
FLAVOR approximation of the solution to the semi-discrete system with the remainders 



(36). 



Proof. The proof of Theorem 5.1 is analogous to that of Theorem 1.2 of [32] (which will 



not be repeated here). The proof requires (46), which is guarantied from Condition 5.2 



by Lemma 5.2 It is easy to check that the slow dependence on time of /, g, rj and fi 
does not affect the proof given in [32]. □ 



Remark 5.3. Condition 1 5 . 2| implies that the constant C in Theorem 5.1 does not depend 
on N or K. This is important because although using a finer mesh leads to a smaller K 
and a larger N = L/K, Condition |5.2| (which is equivalent to the accuracy of the semi- 
discrete approximation of the PDE) ensures that, as long as K = 0{H) and h ^ eH, 



the constant C in the error bounds on the slow component (51 ) and the fast component 



(52) will not blow up. 



Remark 5.4. Observe that the application of the FLAVOR strategy does not require the 
identification of the diffeomorphism 77 (which may depend on the spatial discretization) . 



6 On FLAVORizing characteristics 

The convergence result of the previous section is based on the semi-discretization of 
the original PDE. PDEs and ODEs are also naturally connected via the method of 
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characteristics, and henceforth it is natural to wonder whether a numerical integration 
of those characteristics by FLAVORs would lead to an accurate approximation of the 
solution of the original PDE. The answer to this question will be illustrated by analyzing 
the following (generic) PDE: 

(F{Du,u,q,e-^) = 0, qeU ^^^^ 
\u{q)=-f{q), q£T 

where f7 C M'^ is the domain in which solution is defined, F and 7 define initial/boundary 
conditions. 

The following condition corresponds to assuming that characteristics are well-posed. 
Condition 6.1. Assume that 

1. The PDE F{Du,u,q,e~^) = admits characteristics: 

q = f{q,z,e-') (54) 

z = g{q,z) (55) 

u{q{t)) = z{t) (56) 

where q £ U is a vector corresponding to coordinates of characteristics in the 
domain of the PDE, and z corresponds to the unknown's value along the charac- 
teristics. 

2. For arbitrary e, any point in U is reachable from the initial condition via one and 
only one characteristics. 

The following conditions correspond to the assumption of existence of (possibly hid- 
den) slow and locally ergodic fast variables for those characteristics. 



Condition 6.2. Consider ODE (54). Assume that: 



1. There exists a z-dependent diffeomorphism 77^ : g 1— )• [x,y] from onto "p xW 
with uniformly bounded C^^C"^ derivatives with respect to both q and t, such that 



{x,y) satisfies (with z(t) given by (55 )J 



jx =fi{x,y,z) 

1 • -Iff \ 
[y = e f2{x,y,z) 

where f\, f2, and g have bounded derivatives with respect to x, y and z, and 
u{[r]^]~'^{x,y)) has bounded derivatives with respect to the (slow) variables x 
and z. 

2. There exists a family of probability measures n^{x, dy) on indexed by x £ W^~'p 
and z £ "K, as well as a family of positive functions T 1— )• E^(T) satisfying 
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liiiiT'-^oo E^{T) = 0, such that for all xo,yo, zq,T bounded and (f) uniformly bounded 
and Lipschitz, the solution to 



Yt = f2{xo,Yt,zo) 



In 



yo 



(58) 



satisfies 



(t>(Ys)ds 



(l)iy)n"'ixo,dy) 



xo,yo)\\)E'%T)mL^ + \\V(t>\\L^) 

(59) 

where r i— )■ x^'^ (^) is bounded on compact sets, and /x^ has bounded derivative with 
respect to z in total variation norm. 



The second item of Condition 6.2 corresponds to the assumption that the fast variable 
y is locally ergodic with respect to a family of measures /x drifted by the slow variables 
X and z. 

The following lemma shows that, under the above conditions, the solution of PDE 
(53) is nearly constant on the orbit of the fast components (y) of the characteristics. 



Lemma 6.1. Under Conditions 6.1 and 6.2, for any fixed constant Ci (independent of 
e^"^), there exists a constant C2 independent of e^^ , such that for any < ti < Ci, 
< t2 < Ci and (fixed) xo and zq, 

\u{W']-Hx,,Y{t^)))-u{Wr\x,,Y{t2)))\<C2e 
where Y{ti) andY{t2) are two points on the orbit ofY{t) = /2(xo, l^(t), zq). 



(60) 



Proof. Under Conditions 6.1 and 6.2, it is known (we refer for instance to [29j or to 
Theorem 14, Section 3 of Chapter II of [30] or to [28]) that x and z converge as e — )• 
towards x and z defined as the solution to the following ODEs with initial condition xq 
and Zq 

\^ = J fi{x,y,z)fi^{x,dy) 

= 1 9i[ilT'^{x,y),z)fi^x,dy) 

Therefore, writing y{t) the solution of y = e^^/2(i, y, z), we have as e — >• 

Now, taking the time derivative of n = -u o r/"^, we obtain 

UxX + Uyij + UzZ = z + R{e) 

where R{e) is a function of t that goes to as e — )• 0. 
Furthermore, 



(61) 

(62) 
(63) 



Y{t) = f2{xo,Y{t),ZQ) 

= /2(x(et),y(et),5(et)) + ||(i(et) 
+ o{e) + o{y{et)-Y{t)) 
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By Taylor expansion, x{et) — xq and z{et) — zq are obviously 0{e). Applying Gronwall's 
lemma, we also obtain that y{et) — Y{t) = 0(e). Therefore, 



Yit) = /2(x(et), y(ef), z(et)) + 0(e) = ey(t) + o(e) 
Combining Eq. [63] with Eq. l64l we obtain 



(64) 



t2 



Y{t) dt 



t2 



ydt + o(e) 



UxX — iizz) dt + R{e) + o(e) 



t2 



(65) 



Since u^, x, ut and z are bounded, and R{e) is vanishing (and hence bounded), we 
conclude that the right hand side is 0{e). □ 

Condition 6.3. Assume that the domain U is bounded (independently from e^^). 



Lemma 6.2. // Conditions 6.1, 6.2, and 6.3 hold, then every point in U is reachable by 
a characteristics from the initial condition in bounded time (independently from e~^). 



Proof. From Condition 6.1 , we already know that every point is reachable, and therefore 
it suffices to show that hitting times do not blow up as e — t- 0. Since x{-) converges to 



x{-) (see proof of Lemma 6.1), by considering the x component of the characteristics 



(projected by r/), it becomes trivial to show that the hitting time converges to a fixed 



value (and hence, does not blow up). Using Condition 6.3, we conclude that that any 
point in U can be hit in (uniformly) bounded time from the initial condition. □ 



Analogously to the Integrator 5.1 1 a legacy integrator for (54) and (55) can be FLA- 



VORized, and shown to be convergent under regularity and stability conditions (anal- 
ogous to Condition 5.3) requiring /i, /2 and g to be Lipschitz continuous and q{t) and 



z{t) to be bounded. The convergence result is analogous to Theorem 5.1, modulo the 



following change of notation: the slow index is now z instead of t, the original coor- 
dinates are q instead of Ui, the vector field of the original coordinates is / instead of 
fi, and the dynamics of the slow index comes from the non-trivial drift of i = g{q,z) 
instead of the trivial t = 1. We define u{q{t)) := z{t) for all t on each FLAVORized 
characteristics [q{t),z{t)]. Naturally, u is only defined at discrete points in the domain 
U. These discrete points, however, densely 'fill' the space in the sense that (as shown by 
the proof of the following theorem) FLAVORied characteristics remain very close to ex- 
act characteristics (x components are close in Euclidean distance, and y components are 
close as well in terms of orbital distance induced by the infimum of point-wise Euclidean 
distances). 

By the two-scale convergence theorem, we can quantify: the strong convergence of 
the slow coordinate of the characteristics and the unknown's value along the character- 
istics, and the weak convergence of fast coordinate of the characteristics. Finally, these 
single characteristics' ODE approximation error bounds can be transferred to the PDE 
approximation error bounds by considering the entire family of characteristics starting 
from all points (in initial condition). 
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Theorem 6.1. Write u{q) the solution obtained by FLAVORizing all characteristics. 



Under Conditions 6.1, \6.S\ \6.3[ the consistency and regularity and stability Conditions 



corresponding to Conditions 5.2 and 5.3 (with the change of notation described above), 
there exist a constant C independent of e^^ and qo G T, such that 

\u{q)-u{q)\ <Cxi{qo,j{Qo),e,d,T){l + X2{qo,l{qo),e,S,T-,T,t)) (66) 

for any q on any FLAVORized characteristics, where qo and 7(50) correspond to the 
initial condition that leads to q via a FLAVORized characteristics, and xi o.nd X2 are 
vanishing error bound functions. 

Remark 6.1. When T is compact (such as in the case of periodic boundary condition), 
Xi and X2 can be further chosen to be independent of qo (hence q) by taking a supremum 
over r. 



Proof. By Condition 6.1 all g E U can be traced back to qo & T through a charac- 
teristics. By Lemma 6.2, characteristics starting from qo reach q in bounded time T. 
Using the two-scale convergence of the FLAVORization of these characteristics (a result 
analogous to Theorem |5.1[ ), we deduce that the approximation error associated with zt 
(on FLAVORized characteristics) can be bounded Cxi (with respect to the true value 
u{q) = Zt, the error Ce^^ has been replaced by C because T is bounded). 

Now observe that ijT 7^ qT, where (Jt is the coordinate of the FLAVORized char- 
acteristics starting from qo. As before, let [xT,yT] = vilT) and [xT,yT] = vi'lT)- The 
error on the slow component is — StII < Cxi- The possible large error on the fast 
component is not a problem because we can look for a near-by point on the fast orbit 



with introducing only an 0{e) error on the unknown's value (Lemma |6.l[ ): 

(u{r]{xT, Vt)) = u{r}{xT, y^)) + C'(e) ^^^^ 
\y*T = argmin^^|^^^^(^.^_y^) \\yT - Yt\\ 

Since — a^t|| is small, the local ergodic measures that represent the orbits given by 
Yt = f{xT,Yt) and Yt = f{xT,Yt) will be smaU: \\fi{xT,dy) - /i(xT, dy) ||t.v. < CxiX2 is 
by chain rule. Because yx is on the orbit of Yt = f{xT, Yt), we will have \\y^ — ??^(^t)|| < 
CxiX2- 

All together, we obtain 

lUiqr) - u{qT)\ = \ZT - u{qT)\ 

< \ZT - u{q)\ + \u{qT) - u{qT)\ 

< Cxi + C\\V{u o r?)|U (||XT - V^{qT)\\ + WvT - r?^(9»||) 

< Cxi + C{xi + X1X2) = Cxi + CxiX2 (68) 

□ 

Remark 6.2. To keep the presentation concise, we have written C all constants that 
do not depend on essential parameters. 
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Remark 6.3. As shown above, u will be captured strongly. Du, on the other hand, 
depends on a derivative with respect to the fast variable, and therefore will only be 
convergent in a weak sense. 

Relevance to an error analysis for PDE-FLAVORS The above result guaran- 
tees the convergence of FLAVORized characteristics. It is also possible to establish an 
error bound on the difference between a specific PDE-FLAVOR discretization and the 
approximation given by the above FLAVORized characteristics (and hence prove the 
convergence of this specific PDE-FLAVOR discretization). Such an error bound could 
be obtained by first transforming FLAVORized characteristics to PDE-FLAVOR grid 
points via interpolating functions, and then using the fact that coordinate transforma- 
tions do not affect the efficiency of FLAVORS. For the sake of conciseness, we did not 
elaborate on this point here. 
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